home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / poly / dd.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-12-02  |  2.3 KB  |  102 lines
  1. /* interpolation/interp_poly.c
  2.  * 
  3.  * Copyright (C) 2001 DAN, HO-JIN
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Modified for standalone use in polynomial directory, B.Gough 2001 */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_errno.h>
  24. #include <gsl/gsl_poly.h>
  25.  
  26. int
  27. gsl_poly_dd_init (double dd[], const double xa[], const double ya[],
  28.                   size_t size)
  29. {
  30.   size_t i, j;
  31.  
  32.   /* Newton's divided differences */
  33.  
  34.   dd[0] = ya[0];
  35.  
  36.   for (j = size - 1; j >= 1; j--)
  37.     {
  38.       dd[j] = (ya[j] - ya[j - 1]) / (xa[j] - xa[j - 1]);
  39.     }
  40.  
  41.   for (i = 2; i < size; i++)
  42.     {
  43.       for (j = size - 1; j >= i; j--)
  44.     {
  45.       dd[j] = (dd[j] - dd[j - 1]) / (xa[j] - xa[j - i]);
  46.     }
  47.     }
  48.  
  49.   return GSL_SUCCESS;
  50. }
  51.  
  52. #ifndef HIDE_INLINE_STATIC
  53. double
  54. gsl_poly_dd_eval (const double dd[], const double xa[], const size_t size, const double x)
  55. {
  56.   size_t i;
  57.   double y = dd[size - 1];
  58.  
  59.   for (i = size - 1; i--;)
  60.     {
  61.       y = dd[i] + (x - xa[i]) * y;
  62.     }
  63.  
  64.   return y;
  65. }
  66. #endif
  67.  
  68. int
  69. gsl_poly_dd_taylor (double c[], double xp, 
  70.                     const double dd[], const double xa[], size_t size,
  71.                     double w[])
  72. {
  73.   size_t i, j;
  74.  
  75.   for (i = 0; i < size; i++)
  76.     {
  77.       c[i] = 0.0;
  78.       w[i] = 0.0;
  79.     }
  80.  
  81.   w[size - 1] = 1.0;
  82.  
  83.   c[0] = dd[0];
  84.  
  85.   for (i = size - 1; i > 0 && i--;)
  86.     {
  87.       w[i] = -w[i + 1] * (xa[size - 2 - i] - xp);
  88.  
  89.       for (j = i + 1; j < size - 1; j++)
  90.     {
  91.       w[j] = w[j] - w[j + 1] * (xa[size - 2 - i] - xp);
  92.     }
  93.  
  94.       for (j = i; j < size; j++)
  95.     {
  96.       c[j - i] += w[j] * dd[size - i - 1];
  97.     }
  98.     }
  99.  
  100.   return GSL_SUCCESS;
  101. }  
  102.